#15
X <- scan() # Read the data below in object 'X'
10 15 78 2 98 35
27 22 47 35 102 4
5 7 12 66 190 34
23 17 22 24 56 13
9 12 7 32 41 19
21 17 415 47 70 24
11 4 12 16 23 87
12 67 13 21 8 18
83 38


X_ <- X # copy of X in X_

N <- length(X_)
n <- 10

set.seed(2012)
U <- runif(N)
P <- n*X_/sum(X_)
R <- U/P*(1-P)/(1-U)
which(R<0) # Unit 33 will always be selected


X_ <- X_[-33]
N <- length(X_)
n <- n-1

set.seed(2012)
U <- runif(N)
P <- n*X_/sum(X_)
R <- U/P*(1-P)/(1-U)
which(R<0) # Unit 17 will always be selected


X_ <- X_[-17]
N <- length(X_)
n <- n-1

set.seed(2012)
U <- runif(N)
P <- n*X_/sum(X_)
R <- U/P*(1-P)/(1-U)
which(R<0)
round(R,2)


id <- (1:N)[(R<=sort(R)[n])*(1:N)>0]
id[id>31] <- id[id>31]+1
id[id>16] <- id[id>16]+1
id <- c(17,33,id)

ppips <- rbind(id,c(X[c(17,33)],X_[(R<=sort(R)[n])*(1:N)>0]))
rownames(ppips) <- c("Id","Xi")
ppips




#16a
x16 <- matrix(c(592,22,614,298,88,386,890,110,1000),3,3)
n16 <- x16[3,3]
N16 <- 6696530
p_a <- x16[3,1]/n16
sep_a <- ((1-n16/N16)*p_a*(1-p_a)/(n16-1))^.5

p_a
sep_a



#16b
pc <- .912
p_b <- x16[1,1]/x16[1,3]*pc+x16[2,1]/x16[2,3]*(1-pc)
sep_b <- ( (1-n16/N16)*(pc*x16[1,1]/x16[1,3]*x16[1,2]/x16[1,3]/(x16[1,3]-1) + (1-pc)*x16[2,1]/x16[2,3]*x16[2,2]/x16[2,3]/(x16[2,3]-1)) )^.5

p_b
sep_b



#16c
x16_ <- array(c(355,15,370,155,25,180,510,40,550,237,7,244,143,63,206,380,70,450),c(3,3,2)) 
py <- .421
pc_ <- c(pc,1-pc)
py_ <- c(1-py,py)

#raking
5 -> antal
raking <- array(NA,c(3,3,antal*2+1))
raking[1:2,1:2,1] <- x16_[1:2,3,]
raking[1:2,3,1] <- rowSums(raking[1:2,1:2,1])
raking[3,,1] <- colSums(raking[1:2,,1])

is.even <- function(x){ x %% 2 == 0 }  # create function to check if even or odd

# do the raking
for (i in 1:(2*antal)) {
  if (is.even(i)==F) {
    raking[1:2,1:2,i+1] <- raking[1:2,1:2,i]*pc_*n16/raking[1:2,3,i]
    raking[1:2,3,i+1] <- rowSums(raking[1:2,1:2,i+1])
    raking[3,,i+1] <- colSums(raking[1:2,,i+1])
  }
  if (is.even(i)==T) {
    raking[1:2,1:2,i+1] <- t(t(raking[1:2,1:2,i])*py_*n16/raking[3,1:2,i])
    raking[1:2,3,i+1] <- rowSums(raking[1:2,1:2,i+1])
    raking[3,,i+1] <- colSums(raking[1:2,,i+1])
  }
} # end

#poststratification
post <- rbind(x16_[1:2,1:2,1]/x16_[1:2,3,1]*raking[1:2,1,antal*2+1],
x16_[1:2,1:2,2]/x16_[1:2,3,2]*raking[1:2,2,antal*2+1])

p_c_ <- post[,1]/rowSums(post)
p_c <- sum(p_c_*rowSums(post))/n16
sep_c <- (sum( (rowSums(post)/n16)^2 * p_c_*(1-p_c_)/(rowSums(post)-1)*(1-n16/N16)))^.5

p_c
sep_c